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We discuss an algorithm for the exact sampling of vectors v £ [0,1]^^ satisfying a set of pairwise 
difference inequalities. Applications include the exact sampling of skew Young Tableaux, of con- 
figurations in the Bead Model, and of corrugated surfaces on a graph, that is random landscapes 
in which at each vertex corresponds a local maximum or minimum. As an example, we numeri- 
cally evaluate with high-precision the number of corrugated surfaces on the square lattice. After an 
extrapolation to the thermodynamic limit, controlled by an exact formula, we put into evidence a 
discrepancy with previous numerical results. 



I. INTRODUCTION 

Consider the following problem: given a grid L x L, in 
how many ways can we fill the boxes with the numbers 
from 1 to in such a way that odd/even boxes are 
local maxima/minima? A typical allowed configuration 
for L = 5 is 
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A combinatorial more general setting is the following. 
Consider a direct acyclic graph G = (V, E) , with N ver- 
tices: how many one-to-one maps a : V{G) ^ {1, . . . N} 
exist such that, for each oriented edge {ij) G E{G), 
a{i) > The case of the square grid pertinent to 

Corrugated Surfaces corresponds to the graph 




in ©AT, and the fraction of maps satisfying all the con- 
straints is given by 

- ^ #{^ e &N ■■ m) e E, ail) > aij)} . (1) 

There are many other applications in Combinatorics and 
Statistical mechanics. A connection is with the counting 
of acyclic orientations of a graph, a special evaluation of 
its Tutte polynomial, which is a worst-case #P problem 
(see e.g. S]). 

For a given (unoriented) graph G = (V,E), an acyclic 
orientation (/> is a choice of orientation for the edges, such 
that no oriented cycles are created. Call A{G) the set of 
such 0's, and, for any (j), call G{(j)) the corresponding ori- 
ented acyclic graph. Any non-degenerate function r from 

to a totally ordered set induces an acyclic orientation: 
just orient the edge (ij) from i to j if r(i) > r(j). Call 
G[t] the induced oriented acyclic graph. Given a what- 
ever ordering of the vertices, permutations are special 
cases of valid functions r. 

The uniform measure over ©^r is easy to study an- 
alytically, or to exactly sample. However, it induces a 
biased measure over A{G) (through the natural function 
4>{t) as the orientation (f> such that G[t] = G{(f)). The 
corresponding bias factor is exactly -^g^V]' that, for 
example. 
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This wider framework includes, among other things, the 
counting of Standard Young Tableaux and Standard 
Skew Young Tableaux [U, as in the example below 
with the skew tableau (6,4,4,3) \ (3,2): 
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Of course, given a whatever ordering of the vertices in 
V{G), maps cr are naturally identified with permutations 



and in particular the average of •^g^'r] is related to the 
cardinality of A{G). 

Label the vertices with indices from 1 to A^. An equiv- 
alent formulation of the problem is to ask for the Les- 
besgue measure, in the interval [0, 1]^, of the vectors 
X — {xi\i<i<:M such that Xi > Xj for each oriented edge 
(ij) in the graph, oriented from i to j. Indeed, at the 
aims of A^-dimensional measure, we can neglect config- 
urations with repeated entries, and the constraint only 
depends on the ordering of the variables, so that 
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(here e{x) = 1 if x > and e{x) = if a; < 0). This 
alternate perspective gives justice to the name of "Cor- 
rugated Surfaces" for the configurations on the grid (see 
Figure [1]), and makes expHcit another speciahzation of 
the generahzed model, the "Bead Model" y , correspond- 
ing to a realization on a square lattice rotated by an angle 
7r/4, and the edges being directed, say, in the down-right 
and down-left directions. 

The more general problem of understanding the statis- 
tics of local minima in a random landscape defined on 
a graph, which arises in the study of a point particle 
in a random potential or, for example, by looking at 
the energy landscape corresponding to configurations in 
many-body systems, was addressed in In particu- 

lar in fo'], instead of looking at typical configurations, the 
interest is shifted to the probability of large deviations, 
that is of configurations which are far from typical, as 
the configurations in which the local minima are max- 
imally packed, i.e. Corrugated Surfaces. They provide 
both analytical (for Cayley trees) as numerical estimates 
(for hypercubic lattices) for the constant 7 which con- 
trols the number of those configurations in the limit of 
large iV, as Zq ^ ■ But the numerical results don't 
rely on Monte Carlo simulations. 

In this paper we present a Monte Carlo algorithm to 
"exactly" , that is without any bias, sample configurations 
X G [0,1]^ satisfying our constraints, for an arbitrary 
given oriented acyclic graph G. 

The algorithm runs under the paradigm of Wilson 
and Propp "CoupUng From The Past" (CFTP) 0, in a 
variant which allows for continuous- valued variables (the 
original CFTP setting is described for a discrete config- 
uration space). As an application, we calculate numeri- 
cally the large-size asymptotics (for N — > 00) of the num- 
ber of corrugated surfaces in two dimensions, with four 
digits of precision, under some hypothesis on the scaling 
exponents of the finite-size corrections (that we motivate 
theoretically in the following). 

At the light of this observation, it is not surprising 
that Corrugated Surfaces and the Bead Model allow for 
exact sampling using CFTP, as the relative discretized 




FIG. 1: An example of corrugated surface in a 8 x 8 square 
geometry, and "continuous" formulation as in equation 



variants, hard-core lattice gas on the square lattice and 
lozenge tilings, are among the most studied applications 
of discrete CFTP or similar techniques (see Q for the 
hard-core gas, plus J9|l for a discrete-variable continuous- 
space version, and [lO| for Lozenge tilings). 

The extension of the method to continuous variables is 
not a big deal, and we do not claim much originality in 
this. However it is not either obvious a priori, as a crucial 
point is the possibility that two coupled random Markov 
chains reach coalescence, an event that could naively be 
though to have zero probability in Lesbesgue sense. Some 
of these aspects have been already considered in the liter- 
ature while issues of different sort, and more specific 
to our problem, are discussed in detail in the following. 



II. COUPLING FROM THE PAST 

As we said in the introduction, we used an "exact sam- 
pling" algorithm, i.e. an algorithm that allows to extract, 
in a finite time, a configuration in an ensemble with a 
given measure, without any bias. For a general discus- 
sion on Monte Carlo methods in Statistical Mechanics 
we refer to fl^l where both the initialization bias, that 
is the source of systematic error related to "thermaliza- 
tion" , which is controlled by the exponential autocorrela- 
tion time, and the statistical error, controlled by the inte- 
grated autocorrelation time at "equilibrium" are fully dis- 
cussed. See also [31 for an example of rigourous bounds 
on the autocorrelation times related to physical observ- 
ables for Metropolis algorithms. 

The CoupUng From The Past method (CFTP) also 
makes use of Markov Chain processes. However, the two 
fundamental concepts of evolving a set of coupled chains, 
and analysing the dynamics backward in time, allow to 
determine a paradigm which leads to perfectly unbiased 
samples 0. 

A coupled Markov chain is a process involving a set 
of configurations (say, k of them), which evolve under 
a dynamics being a valid (ergodic equilibrium) dynam- 
ics for each chain, but using the very same sequence 
of randomly-generated numbers in the different copies. 
This makes it a non-ergodic dynamics for the A:-uple as 
a whole. In particular, if any two configurations become 
the same at some time t, they will not anymore evolve 
into distinct configurations. The first time at which all k 
chains reach the same configuration is called coalescence 
time, and the coupled chain is said to have coalesced in 
this case. 

Consider a conceptual experiment in which we follow 
the evolution of all the allowed configurations, with the 
chains weighted according to the measure on the ensem- 
ble. We have that, at least at the initial time, averages on 
the k copies of the chain correspond to statistical averages 
in the ensemble. Furthermore, this property must be pre- 
served by the time evolution, and in particular must be 
true also in a limit i — > 00, in which, if the average time 
of coalescence is finite, we have a single configuration. 
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In order not to introduce a bias due to the "obser- 
vation" that the chains have coalesced, one must use a 
protocol of "looking backward" in time, i.e. imagine to 
perform the coupled-chain experiment above starting at 
a sequence of negative times, for example tg = — 2^, and 
stopping the dynamics at t — 0. Each finite state of the 
coupled chain has the property that averaging on all the 
survived different copies, with the measure induced by 
the (averaged) time evolution, would provide the appro- 
priate statistical averages, and in particular, at the first 
generation g such that the chain has coalesced, the only 
survived final state is exactly sampled. 

Following the full set of possible states is clearly un- 
feasible. However, there is a case in which one can cer- 
tificate that all the states in the conceptual experiment 
above have coalesced, by just following a coupled Markov 
Chain with k of order 1. This happens if the space of 
configurations has a structure of (mathematical) lattice, 
i.e. there is a structure of partial ordering ^, and two 
special states O and / such that, for each X in the lat- 
tice, O < X ^ I Furthermore, it is required that 
the (ergodic reversible) Markov Chain dynamics, whose 
equilibrium distribution is the desired measure on the en- 
semble, preserves the ordering, i.e., the coupled evolution 
of two ordered configurations Xt d: Yt leads to ordered 
configurations Xt+i ^ Yt+i. 

In this case, it suffices to work with k — 2, and initial 
states O and /, in order to have that if O and / have 
coalesced, also all configurations in between did it. 

In summary, we can extract a general "CFTP pro- 
tocol" for sampling on a lattice. Given an ensemble, 
suppose that you can find a lattice structure on the 
space of configurations, and an ergodic reversible dy- 
namics which preserves the ordering. Then, iteratively 
for (7 = 0, 1,2,..., run the k = 2 coupled Markov Chain 
with initial states {O, I). It is important that, in the last 
2^ times of the simulation at generation 5 -I- 1 you use 
the same set of random numbers used in all the 2^ times 
of the run at generation g. Stop the simulation at the 
first generation g^, such that the chain has reached coa- 
lescence: the exactly sampled configuration is the state 
at time 2^* . 



III. A CFTP ALGORITHM FOR OUR MODEL 

Given our directed acyclic graph G — (V^E), for each 
vertex i we define J^+{i) as the set of vertices j such 
that (ji) G E{G), and J\f-{i) as the set of vertices j such 
that {ij) G E{G). Then, a restatement of the constraint 
n(ij) ^(-^i ~ -^j) that, for each vertex i e V{G), 



max (xi) < Xj < 



min (xn 



(4) 



It is easy to devise an ergodic Monte Carlo chain with 
single- variable heat-bath moves. More explicitly, at each 
time step 



1. Choose i at random uniformly in V{G); 

2. Choose z at random uniformly in [0, 1]; 

3. Replace Xi by z if it happens that 

max (xi) < z < min (xj) . 

Remark how conversely, in the discrete formulation, such 
a dynamic would have been "quenched" by the further, 
highly non-local constraint that all (t(i)'s are distinct. 

The simple Monte Carlo chain above has the remark- 
able property of being suitable for CFTP, in the way de- 
scribed in Section ini this being another advantage of the 
continuous formulation, w.r.t. the one in terms of permu- 
tations. Consider the space S C [0, 1]^ of valid vectors 
X, under the natural partial ordering x ^ y ii xi < yi for 
all i G V{G). For any graph G, this space is a Lattice, 
as there are both a O and a / element, corresponding to 
the vectors and 1 respectively. 

Then, the second condition for CFTP is that, given 
two configurations x and y at time t such that x the 
coupled time evolution preserves the ordering. This is 
easily seen, with the help of the crucial observation that, 
for any vertex z, with the definitions 
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argmax(xj ) ; 
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one has 



min {x.j) < Xj. < yf = min {y,) ; 
max {xj) = Xj>> < yj» < max (y^) 

jGA'-(i) je7V_(i) 



Call 



x+ = min (xi) ; 



max (xi) ; 



(5) 

(6a) 
(6b) 

(7) 



The range of values in which the new candidate variable 
z is accepted in Xi is the interval and similarly 

for y. The statement in ^ is that x± < y± li x < y. 
If x+ > y-, there is a probability {x^ — y-)/[y+ — x-) 
that, given that a move occurs, the number of indices i 
for which Xi = yi increases by one. If instead Xi = yi at 
some time (and in this case it must he x+ >?/_), there is 
a probability 1 — {x+ — y-)/(y+ — a;_) that this number 
decreases by one. Calling C{x,y) — {i : Xi = yi}, we can 
define a distance parameter d{x,y) = \C{x,y)\, and the 
reasonings above describe a non-trivial hopping dynamics 
in the parameter c?(x, y) S {0, . . . , N}. The chain starts 
from d{x, y) — N and reaches coalescence when d{x, y) = 
0, which is thus a fixed point of the induced restricted 
dynamics. 

One could naively think that the extrema d = and 
d = N are strongly repulsive, in a way that becomes 
stronger with lattice size, because of "entropic reasons" 
(e.g., for d = 0, if |C| = 0{N) and \V{G) \ C| = 0(1), 
we have a relative factor 1/iV for choosing an index in 
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the second set). At least in graphs with low degree, say 
bounded by k, this is the case for the d — N, but not 
for d — 0. Indeed, if the dynamics chooses a site in C 
such that all its neighbours are in C, the distance can 
not increase at that time step, so that the ratio between 
the two rates d d ±1 is bounded by k. The underly- 
ing mechanism is somehow related to the fact that naive 
entropic reasonings hold for equilibrium statistical me- 
chanics, but require more care in non-ergodic dynamics 
such as in the coupled Markov Chain. 



IV. FROM EXACT SAMPLING TO THE 
CALCULATION OF THE FREE ENERGY 

It is well known [3] (see also [iB, sec. 3.2]) that the prob- 
lems of (approximated) counting and of uniform sampling 
are closely related. In particular, there is a way of calcu- 
lating numerically the free energy of a model of interest 
if we can generalize the model introducing an extra pa- 
rameter e such that 

• at e = 1 we recover the model of interest; 

• at e = the free energy is known exactly; 

• exact sampling is available in the range e e [0, 1]; 

• the (unnormalized) Gibbs measures /i^, for differ- 
ent values of e, are absolutely continuous (in either 
direction). 

Actually, instead of the last point, it is sufficient to have 
the weaker but more technical statement of having a 
Radon-Nikodym derivative of a measure at a value e 
w.r.t. a measure at a value e', for a suitable set of pairs 
(e, e') (see for the pertinent definitions). Indeed, if, 
for example, for e < e' we have /i^ <IC fi^' , then we have a 
fimction g^^^i{x) such that 



(8) 



and thus, by defining the free energy as F{e) = In | (we 
neglect overall constants customary in thermodynamics) , 
we have 



F(e)=i^(e')+ln(5e,e').. , 



(9) 



which can be used telescopically in order to obtain F{1) 
from F{0) and the evaluation of In (g^.^e.^^) through the 
exact-sampling algorithm run at some sequence of Ci's, 
sufficiently dense that the statistics of g^^^^i-^i is significa- 
tive. Also, if one can define 



lim 

5^0 



In 



(10) 



and this function is smooth in e, the free energy for the 
model of interest would be given by 



and it could be more efficient to find a reasonable con- 
tinuous fit of g{e) from simulations run at some sequence 
of Ei's. 

In our case of corrugated surfaces, for a configuration x 
define W[x\ C V{G) as the set of vertices being maxima 
and in the range [0, union the ones being minima and 
in the range Call m{x) = \ W[x\\. Then choose 



IJ.,{x) 



m{x) 



[12) 



It is evident that t — 1 corresponds to our model, and 
that at e = we just have Zq — 2^^ . 

The Markov Chain introduced in Section IIIII is easily 
generalized to the introduction of the parameter e. Just 
the second point, where we ask to extract z uniformly in 
the interval [0, 1], has to be replaced by the measure on 
[0,1] 

p^^'(^) = YT7G + ('"^)'(^(^"^))) 

with +/— if we are performing the move respectively on 
the position of a maximum or a minimum. The plot of 

p^e^\z) is just 
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Specialization of the quantity in (|10p is easily achieved: 



9e 



•(^) - (?) 



and, taking the limit. 



g{() = - {m{x))^ 



(13) 



(14) 



The apparent singularity at e = is not there, as (m{x))^ 
vanishes linearly with e. Also, the potential risk of having 
an explosion of data noise at e = is easily avoided: as 
this point is the trivial limit of the theory, it is easy to 
match the data near e = with the first few terms of a 
Cluster Expansion. Call G{L) the graph corresponding 
to a grid of side L, and define the intensive free energy 
as 
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InZ, 



G(L) 



(15) 



then the Cluster Expansion, in powers of e and inverse 
powers of L, gives 



F{1)^F{0) + J\eg{e), (H) = " 2 + ^ - (^^) + ^, ^) . (16) 
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FIG. 2: Distribution of L'^m{x), averaged over 10'' exactly 
sampled configurations, on a system with size L = 256, at 
e — 1. The fit is with a Gaussian. 



The observable W[x] has also an appealing interpreta- 
tion. It is clear that no adjacent vertices can be simulta- 
neously in W, so, for each configuration x, W is a hard- 
core gas configuration (or an independent set). The pa- 
rameter e plays the role of an effective fugacity in the 
gas, and, although the correspondence is not perfect, we 
can imagine that, if any criticality at all appears in the 
model, it will be in the same universality class of the one 
of the two-dimensional hard-core lattice gas. 

For the general case, one can alternatively use the intu- 
itive "temperature" parametrization. Allow for all con- 
figurations X G [0, 1]''^, and define h{x) as the number 
of edges whose constraint is not satisfied. Then we can 
write 



fXeix) 



n 

(y)6_E 



eeix. 



(17) 



Again e = 1 corresponds to our model, while at e = we 
just have Zq = 1. Similarly, 



9e,e'{x) 



and, taking the limit. 
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1 - e' 



h{x) 



m^-j—^{h{x))^ 



(18) 



(19) 



Again the apparent singularity at e = 1 is not there, 
as (h{x))^ vanishes linearly at that point. However, in 
this case the increase of noise at e = 1 is unavoidable, 
this being the reason why we have chosen the ad hoc 
parametrization ()12[) for our numerical simulations of cor- 
rugated surfaces. 



V. NUMERICAL DATA FOR CORRUGATED 
SURFACES 

In order to obtain a reliable estimate we performed our 
simulations for values of L = 16, 32, . . . , 256, and at 
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FIG. 3; Numerical data for the function in ([22}. The natural 
error bars are invisible. In the figure, we magnified the errors 
by a factor 100 in order to highlight the different relative 
error in the two regimes of e — > and e = 0(1). The linear 
behaviour at e = is deduced from the Cluster Expansion 
in (fT6l). 



20 different values of e, equally spaced, in the relevant 
interval, in the parametrization e/(l -t- e). For each value 
of L and e, we performed 10* independent runs, providing 
us a set of numerical values for the quantities 



(20) 



For each value of e in our analysis, the distribution of 
m(x) in the output configurations is well-fitted by a 
Gaussian (cfr. for example Figure [5]), in agreement with 
the fact that, if rn(a;) is a good order parameter, as we 
expect from the analogy with the hard-core lattice gas, 
we are in a regime with a single Gibbs phase. We checked 
that our numerical implementation provides exactly the 
results on the non-trivial case L ~ 3 which can be com- 
puted analytically. For the larger values of L, at each 
value of e, the numerical values of ^/(e; L) are in good 
agreement with a finite-size description in which the first 
correction scales with 1/L. So, in a fit of the form 
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we extrapolated the asymptotic values 
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(21) 



(22) 



A scrupolous statistical analysis of the errors, just within 
the (well-verified) assumption on the exponent of the 
leading finite-size correction, leads to the results in Fig- 
ure [3l Finally, numerical integration of a polynomial in- 
terpolation of the data (with a polynomial of degree 5, 
determined by analysis of structure in the errors) pro- 
vided us with the result 
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-0.53967 ±0.00054. 



(23) 
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The parameter 7 is obtained from 7 = exp(— /). Our 
numerical analysis on the square lattice gives 



g : 14 15 16 17 



7 = 1.7154 ±0.0009 



(24) 



with purely statistical errors. This must be compared 
with the previous numerical estimate ^ 



7mm(2) = 1.6577 ±0.0006, 



(25) 



which is definitely out of the estimated errors. Even if the 
difference appears to be small it is interesting to observe 
that, as the estimate in 6] for the three-dimensional cubic 
lattice is 



7mm(3) = 1.7152 ±0.0010 



(26) 



the discrepancy is of the same size of the estimated differ- 
ence between the two and three dimensional constants. 



VI. ANALYSIS OF THE TIMES TO 
COALESCENCE 

According to the "ordinary" CFTP protocol (not the 
Read-Once protocol of [§]), the running times for each 
size, and each independent instance, is essentially a power 
of 2, say, 2^. This means that the average complexity of 
the algorithm is described through some probability for 
the generation parameter g, at side L, that we call Phig) 



Time(L)^^29pL(g) 



(27) 



It is a tautology that this time is finite provided that 
Phig) has a Laplace Transform at the value — In 2. It 
comes out that, in our range e € [0, 1], the times grow 
monotonically with e, and even at the hardest point e = 1 
the histogram Phig) takes its leading contribution from 
one or two values of 5. At fixed g, the plot of pl (g) in the 
variable L instead looks like a smooth bell-shaped curve 
in the range [0, 1] (cfr. Figure [3]), and it is fairly safe to 
interpolate the (continuous) value L{g) of L at which the 
curve for g and the one for g + 1 do cross. Given this 
regular behaviour, the scaling of L{g) with g must be 
related to the scaling of the complexity, in particular, if 
g{L) is the functional inverse of the (obviously monotone) 
L{g), then 

Time(L) - 2^^^^ . (28) 

The data are well fitted by a curve of the form 

g{L) (xL'^ In L, (29) 

so the complexity seems to be only logarithmically super- 
linear in the number of degrees of freedom. 

As we discussed in Section ITVl the pictorial interpre- 
tation of the observable W[x] suggests that our model of 
Corrugated Surfaces, in its continuation to arbitrary val- 
ues of e, is in the same universality class as the hard-core 
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FIG. 



4: On the left, plots of pL{g) in the variable L, for g = 
8, . . . , 19. On the right, table of the values L{g) interpolated 
from the data. Errors are estimated to be on the last digit. 



lattice gas, so it was natural to expect that, analogously 
to what is proven for this model [l^, [l^, the average 
coalescence time is ~ n In n for Glauber dynamics at suf- 
ficiently small fugacities (in agreement with ((29)l above), 
while it becomes worst-case hard at sufficiently large fu- 
gacities, these values being lower- and upper-bounds to 
the "physical" critical values e* (the gap being originated 
by technicalities in the proof procedure) . It is not incon- 
ceivable that, for our system, all the range e < 1 has fast 
coalescence times. 

Furthermore, fast convergence would imply a choice 
of parameters which are far from a critical point, and 
correspond to a single thermodynamic phase [20, 21], and 
this would imply in turns that finite-size corrections to 
intensive observables scale with the ratio perimeter/area, 
i.e. with L^^ in our two-dimensional case. This justifies 
the treatment of the finite-size corrections that we have 
done in Section fVl 



VII. CONCLUSIONS 

We have pointed out how the method of Coupling From 
The Past may be fruitfully applied to problems with vari- 
ables assuming values on a continuum domain, and how 
a one-variable Heat-bath Monte Carlo chain is suitable 
for CFTP in an interesting problem in this class: the 
uniform sampling of height functions satisfying a set of 
inequalities described by an acyclic graph. 

As a specific example, we performed a numerical sim- 
ulation for the model of Corrugated Surfaces, in order 
to determine the value of the free energy. Our preci- 
sion goal was to have error bars much below the order 
of magnitude of the discrepancy between our and the es- 
timate, obtained without using a Monte_ Carlo method, 
that appears in Majumdar and Martin [6]. To this aim, 
the use of an exact sampling method in connection with 
an exact formula for the extrapolation to the thermody- 
namical limit has been useful in order to rule out any 
possible source of systematic errors. 

We have achieved our goals with a relatively small nu- 
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merical effort, because the model of Corrugated Surfaces 
appears to be a special value of a one-parameter family of 
models, in the universality class of the two-dimensional 
hard-core lattice gas, in the low-density phase. In this 
case, as expected from the literature, the mixing times 
are only logarithmically super-linear, and the coalescence 
time for the couple chain is of the same order of magni- 
tude of the mixing time of the ordinary Markov Chain. 
As a final comment, we advice that, in this regime, the 
strong control over the errors in the CFTP overwhelms 



the small gain in terms of computational times of other 
Monte Carlo algorithms. 
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